# -*- coding: utf-8 -*-
"""
    @project: pythonProject
    @Author：HanYonghua
    @file： fir1.py
    @date：2025/7/4 12:54
    @blogs: https://www.ncatest.com.cn
"""

import numpy as np
from scipy.signal import freqz
import matplotlib.pyplot as plt


def design_fir_lowpass(cutoff_hz, fs, num_taps, window='hamming', scale=1.0):
    nyquist = fs / 2
    normalized_cutoff = cutoff_hz / nyquist
    n = np.arange(num_taps)

    # 可调缩放系数
    h_ideal = 2 * normalized_cutoff * np.sinc(scale * normalized_cutoff * (n - (num_taps - 1) / 2))

    if window == 'hamming':
        win = np.hamming(num_taps)
    elif window == 'blackman':
        win = np.blackman(num_taps)
    else:
        win = np.ones(num_taps)

    h = h_ideal * win
    h = h / np.sum(h)
    return h

def design_fir_lowpass2(cutoff, num_taps, window='hamming', scale=1.0):
    normalized_cutoff = cutoff
    n = np.arange(num_taps)

    # 可调缩放系数
    h_ideal = 2 * normalized_cutoff * np.sinc(scale * normalized_cutoff * (n - (num_taps - 1) / 2))

    if window == 'hamming':
        win = np.hamming(num_taps)
    elif window == 'blackman':
        win = np.blackman(num_taps)
    else:
        win = np.ones(num_taps)

    h = h_ideal * win
    h = h / np.sum(h)
    return h

# 参数设置
fs = 1000
cutoff_hz = 100
num_taps = 201
window = 'hamming'
scale = 1.02  # 微调系数

# 设计滤波器
h = design_fir_lowpass(cutoff_hz, fs, num_taps, window, scale)
h = design_fir_lowpass2(0.2, num_taps, window, scale)



# 验证
w, H = freqz(h, 1, worN=8000, fs=fs)
H_db = 20 * np.log10(np.abs(H))
cutoff_idx = np.argmin(np.abs(H_db - (-3)))
cutoff_actual = w[cutoff_idx]
print(f"实际-3 dB截止频率: {cutoff_actual:.1f} Hz")

# 绘图
plt.figure(figsize=(10, 4))
plt.plot(w, H_db, 'b', label='Frequency Response')
plt.axvline(cutoff_hz, color='r', linestyle='--', label=f'Target: {cutoff_hz} Hz')
plt.axvline(cutoff_actual, color='g', linestyle=':', label=f'Actual: {cutoff_actual:.1f} Hz')
plt.xlim(0, fs / 2)
plt.ylim(-120, 5)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Gain (dB)')
plt.title(f'Optimized LPF: {cutoff_hz}Hz, Scale={scale}')
plt.grid(True)
plt.legend()
plt.show()